Possible Exotic phases in the One— Dimensional Extended Hubbard Model 
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We investigate numerically the ground state phase diagram of the one-dimensional extended Hub- 
bard model, including an on-site interaction U and a nearest-neighbor interaction V . We focus on 
the ground state phases of the model in the V S> U region, where previous studies have suggested 
the possibility of dominant superconducting pairing fluctuations before the system phase separates 
at a critical value V — Vps- Using quantum Monte Carlo methods on lattices much larger than 
in previous Lanczos diagonalization studies, we determine the boundary of phase separation, the 
Luttinger Liquid correlation exponent K p , and other correlation functions in this region. We find 
that phase separation occurs for V significantly smaller than previously reported. In addition, for 
negative U, we find that a uniform state re-enters from phase separation as the electron density is 
increased towards half filling. For V < Vps, our results show that superconducting fluctuations are 
not dominant. The system behaves asymptotically as a Luttinger Liquid with K p < 1, but we also 
find strong low-energy (but gapped) charge-density fluctuations at a momentum not expected for a 
standard Luttinger Liquid. 



PACS numbers: 71.10, 74.20.Mn 
I. INTRODUCTION 

One-dimensional (ID) Hubbard models have been 
used to model many quasi-lD systems, n including 
conducting polymers such as polyacetylene,EI and or- 
ganic charge-transfer materials such as TTF-TCNQ or 
(TMTSF) 2 PF 6 E In the simplest form, with only an on 
site interaction U, the ID Hubbard n model has been 
solved exactly using the Bethe AnsatzB However, longer 
range interactions are needed to describe many effects 
observed in real materials, e.g., the excitons found in con- 
ducting polymersJa The conventional extended Hubbard 
model (EHM), which in addition to U includes a nearest- 
neighbor interaction of strength V, is the simplest exten- 
sion that takes into account some of the longer-range 
interaction effects. The EHM Hamiltonian is 

itj i 

+V/; n i n i+l> (!) 

i 

where c\ a creates an electron of spin a on site i, and m a = 

c ia c i<? ( n i ~ n iT+ n «l) are electron number operators. An 
implicit parameter is the filling factor n — N e /N, where 
N e and N are respectively the number of electrons and 
lattice sites. We will henceforth give energies in units of 
the hopping t. 

As one of the basic many-body Hamiltonians in ID, 
the EHM—has been the subject of a large number of 
studies.tjEj Nonetheless, as we shall discuss below, there 
remain important open questions related to the phase di- 
agram at intermediate and strong coupling, where both 
analytical and numerical methods are difficult to apply 
reliably. One of the principal reasons for the existence 
of these open questions is the variety of potential bro- 
ken symmetry fluctuations that can occur in the EHM. 



As the parameters U, V, and n are varied, several dif- 
ferent types of correlations dominate the ground state, 
including spin or charge density wave (SDW/CDW), and 
singlet or triplet superconducting (SC) fluctuations. Of 
course, in a strictly ID model, long-range order that 
breaks a continuous symmetry (i.e., SC or SDW) is not 
possible; however long-range CDW order can occur at 
zero temperature. Further, in some parameter regions, 
the EHM is unstable towards phase separation (PS), with 
the nature of the coexisting phases depending on the pa- 
rameters. 

For small values of the interaction parameters, the low- 
energy excitations of the EHM. can be mapped directly 
to a "Luttinger liquid" (LL)Jlil the unifying framework 
for ID interacting fermion systems obtained from weak- 
coupling renormalization group studies, bosonization, 
and conformal field theory. The general (g-component) 
LL contains q gapless degrees of freedom, and the forms 
of the correlation functions depend on only two param- 
eters for each gapless mode a: A renormalized velocity 
v a and an interaction parameter J^__Epr instance, the 
standard Tomonaga-Luttinger modeHE3 has two gapless 
degrees of freedom, charge and spin, and is thus a two- 
component LL with interaction parameters K p and K a \ 
a similar one dimensional exactly solvable model, due to 
Luther and Emery,L3 has gapless charge excitations, but 
a spin gap, and thus behaves as a Luttinger liquid only in 
the charge sector. The identification of a given model as 
a Luttinger liquid enables (in principle) a straightforward 
numerical investigation of the ground state phases of the 
model. In the case of the integrable standard Hubbard 
model, the Bethe Ansatz equations have been used this 
way to calculate the ground state parameters and Kp for 
all values of the repulsion U and the band fillings For 
more general models, away from weak coupling, there is 
unfortunately no reliable analytic method to determine 
the parameters, although they may in principle be calcu- 
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lated numerically. Assuming the model remains a LL also 
away from the weak-coupling regime, numerical methods 
can be used to calculate the exponents K a , which will 
then give the asymptotic form of the correlation func- 
tions. Importantly, for such calculations to be meaning- 
ful, finite-size and other systematic errors must be care- 
fully analyzed, and in particular, one must ascertain that 
the model is not phase separated in the thermodynamic 
limit. As we demonstrate below, calculating Luttinger 
liquid parameters from finite-size systems can sometimes 
be highly problematic, particularly close to phase sepa- 
ration. 

The existence of several distinct regions and tjagpjpf 
phase separation in the phase diagram of tl^r^MBLJtJ'ti!l 
and other ID strongly-correlated modelsEjt^l is well- 
established. While phase separation typically occurs out- 
side the parameter regions thought to be relevant for 
modeling real materials, for purposes of understanding 
the behavior of any given theory it is imperative that 
all phase-separated regions be identified and, if possible, 
understood prior to carrying out other studies, such as 
determinations of Luttinger liquid correlation exponents. 
In small systems, the signals for phase separation can 
be unclear or ambiguous. For instance, the correlation 
functions calculated for small systems can be mislead- 
ing in phase-separated regions, since the boundary be- 
tween the two phases can be large compared to the sys- 
tem size. This leads to a mixing of correlations from the 
two different phases, which may have quite different prop- 
erties (for example, spin-gapped and non-spin-gapped). 
Nonetheless, in a small system, these correlations may 
not show signs of phase separation until well inside the 
phase-separated region. Sometimes, it may then appear 
that the exponent K p > 1 (which would signal dominant 
superconducting correlations), when in fact the system 
is phase separated in the thermodynamic limit and is 
thus not a uniform Luttinger liquid. Accordingly, one 
of our principal goals in this study has been to under- 
stand thoroughly the phase-separated regions before cal- 
culating Luttinger liquid exponents or other correlation 
functions. 

Recently, the extended Hubbard model in the region 
of parameters V ^> U ~ t has attracted considerable 
interest because of the possibility of a novel supercQiu 
ducting ground state away from half- filing (n < 1) .ETeI 
Although the interactions appear to be purely repul- 
sive when V > U > (we also consider V 3> U with 
U < 0), the gradient of the potential is positive, and 
there is an attractive force between electrons. Pairs might 
then form and could in principle lead to a ground state 
with dominant superconducting fluctuations. Based on 
Lanczos exact diagonalization (ED) results at quarter- 
filling (n — 1/2), it was argued that superconducting 
fluctuations indeed dominate for a substantial range of V 
values, before the system phase separates at a very large 
V = Vps -H'Q Another ED study found similar behavior at 
filling n = 2/3, but no attempt was made to determinft 
whether the system is phase separated at this filling.El 



These previous studies also addressed the question of the 
location of the spin gap boundary. Determining the loca- 
tion of the spin-gapped Luther-Emery phase is important 
to studying the possibility of superconductivity in this re- 
gion, as the presence of a spin gap would be consistent 
with the expected dominant singlet pairing correlations 
when V is large and positive. Superconducting corre- 
lations involving triplet pairs have also been proposed 
as a scenario to explain ED results E Triplet supercon- 
ductinffi.|Correlations are dominant in some regions with 

V < OBEij), but explaining their origin for V \U\ is 
problematic. 

In this paper, we investigate the related questions of 
novel superconducting fluctuations, calculations of K pi 
and the boundaries of phase separation and spin-gapped 
regions using quantum Monte Carlo (QMC) simulations 
of relatively large systems (up to 128 sites). We focus on 
the ill-understood parameter region V 3> \U\ for a wide 
range of fillings. Our results show that phase separa- 
tion extends to much lower values of V than previously 
reported.Lrlj We also find that for negative U the high- 
density phase is not the naively expected one with dou- 
bly occupied sites separated by one site (corresponding 
to half filling in the high-density phase) but is at a lower 
density. A uniform state (which has strong CDW fluctu- 
ations) re-stabilizes as half filling is approached. In most 
of the parameter space we can conclude that K p < 1 for 

V < Vps- In some cases, our results are on the border- 
line K p > 1 at our calculated phase separation boundary, 
but in no case do we find a definite K p > 1 before phase 
separation. We therefore believe that superconducting 
correlations never dominate in the V ^> \U\ region. In- 
stead, for a range of parameters, we find strong charge- 
density fluctuations at a momentum 2/cp < <7 < 4fcp, 
which are not expected in a standard Luttinger liquid. 
Our analysis of the temperature and size dependence of 
these fluctuations shows that they do not correspond to 
gapless excitations, and hence the model remains a Lut- 
tinger liquid in the asymptotic low-energy sector. Never- 
theless, these strong anomalous fluctuations demonstrate 
the appearance of non-Luttinger liquid features in the ex- 
citation spectrum. 

To present our results, we begin in Sec. II by dis- 
cussing briefly the numerical methods we have used to 
study the extended Hubbard model. We point out advan- 
tages and disadvantages of several different QMC meth- 
ods and stress the necessity of comparing and contrast- 
ing their predictions to obtain definitive conclusions. In 
Sec. Ill, we discuss methods to detect phase separation 
in numerical data, emphasizing a number of often over- 
looked subtleties, and present our results for the phase 
separation boundaries of the extended Hubbard model. 
We address the determination of Luttinger liquid expo- 
nents and other correlation functions in Sec. IV. Sec. V 
concludes with a summary of our understanding of the 
"exotic" phases of the extended Hubbard model in the 
region V S> \U\. 
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II. NUMERICAL METHODS 

Using the Luttinger liquid formalism, extracting the 
dominant correlations of ID electron systems is easy in 
principle. Lanczos exact diagonalization of small systems 
can give results for certain observables (e.g., velocities 
and stiffness constants) that can be used to determine 
the asymptotic form of the correlations functions which 
are believed to be less affected by finite-size effects than 
the correlation functions themselves. This has led to an 
upsurge of ED studies of various ID models, including 
many focusing on possibilities of dominant SjUntpjimdaicU 
ing fluctuations close to phase separation.^ qEjO'EHTlHI 
However, as we argue below, the finite-size effects may 
in fact be unexpectedly large in regions near phase sep- 
aration, ft is therefore important to confirm ED results 
using methods that can treat larger system sizes. For 
this purpose, we have used three different QMC meth- 
ods. Here we summarize briefly the basic ideas of these 
techniques, referring readers to the literature for addi- 
tional details. Our main purpose is to make some ob- 
servations concerning the advantages and disadvantages 
of the different QMC methods in studies of the EHM in 
the difficult parameter regime V ^S> \U\. We believe that 
most of this discussion will apply to the region of very 
strong interactions in other models as well. 

The first method is based on the "stochas4ji|G-series ex- 
pansion" (SSE) of the density matrix e ~^?B3^C3 This is a 
generalization of Handscomb's method,t3 applicable for 
a much larger class of lattice Hamiltonians. There are no 
approximations such as Trotter discretization of imagi- 
nary time, but in order to obtain ground state results 
one has to ensure that a large enough inverse tempera- 
ture (3 is used. To check for finite-temperature effects, 
we have carried out calculations for several values of (3. 
In general, we find that (3 ~ 2N — AN is sufficient to give 
accurate values for the ground state parameters for N up 
to 128. 

The second technique is a recently developed variant 
of the continuous-time worldlinc algorithm proposed by 
Prokofev et aLe3 Our version of this method^ uses an 
updating scheme adapted from the SSE algorithm Lthc 
two methods are, in fact, closely related to each otherEH). 
The algorithm is based on the finite-temperature per- 
turbation expansion in the interaction representation, 
around the atomic-limit system with no kinetic energy 
term. This expansion converges for a finite system at 
finite [3. The terms (which are paths in continuous imag- 
inary time) can be stochastically sampled in much the 
same way as is done in the SSE method. We will refer 
to this QMC technique as the interaction representation 
(1R) method, ft is also free from systematic errors. 

Results of SSE and IR simulations are in general in 
complete agreement with each other. Both methods 
operate in the real-space occupation number basis and 
hence suffer from well known "sticking problems" (in- 
ability of the local Monte Carlo updates to evolve the 



real space configurations through states with high poten- 
tial energy) when the interactions are very strong. For 
the EHM, V ~ 10 seems to be the highest accessible 
V in the interesting filling regions (U represents a lesser 
problem, since we consider V 3> U). We have found 
that the autocorrelation times for spin and density cor- 
relation functions at strong interactions are shorter for 
the IR method, in particular close to half filling. All 
results presented here for correlation functions and sus- 
ceptibilities were therefore obtained with that method. 
However, for the ground state energy (i.e., the internal 
energy at sufficiently high (3), the SSE method gives re- 
sults with statistical errors approximately an order of 
magnitude smaller than the IR method (in cases where 
the sticking problems do not make simulations practi- 
cally impossible). This result probably arises because 
the total energy estimator of the SSE method is directly 
related to the stochastically averaged order of the terms 
of the expansion of e~@ H . In the IR method (and other 
worldlinc methods), the energy estimator consists of sep- 
arate diagonal and off-diagonal parts, i.e., the potential 
and kinetic energies are not treated on an equal footing 
(in particular, the diagonal part is typically not spin- 
rotationally invariant). We note that the SSE method 
also has proven superior in energy calculations for spin, 
models, such as the two-dimensional Heisenberg m©del,E3 
for which otherwise very efficient loop algorithmsc^l have 
not given nearly as accurate results. 

The third QMC method we employ is the recently 
developed Constrained Path Monte Carlo (CPMC) al- 
gorithm. This is a Slater-determinant based projec- 
tor method that handles the interactions via a Trotter 
decomposition and a Hubbard-Stratonovich transforma- 
tion, leading to auxiliary fields that are sampled.EJ : EJ A 
constraining wavefunction is used to prevent the fcrmion 
sign problem. In ID, the constraint becomes exact (i.e., 
there are no sign problems). There is a small systematic 
error originating from the Trotter decomposition, which 
can be made arbitrarily small by decreasing the "imag- 
inary time slice" width. In our comparisons with SSE 
and IR results, CPMC yielded similar results, except at 
very strong interactions where matrix conditioning prob- 
lems become overwhelming and make the method very 
hard to use (much before the sticking problems become 
problematic in SSE and IR simulations). One advan- 
tage of CPMC over the SSE and IR methods is that in 
CPMC any equal-time correlation function may be com- 
puted since the Green's function is directly accessable 
in CPMC. To obtain accurate energies from CPMC it 
was necessary to compute the energy for more than one 
At value and scale the results to At — ► 0. Our imple- 
mentation of CPMC used a uniform free-electron wave- 
function for the constraint, the initial wavefunction, and 
the importance sampling wavefunction. While the choice 
of the importance sampling wavefunction should not af- 
fect the final results, if the overall symmetry is incor- 
rect the method becomes very inefficient. In phase sep- 
arated regions, we found that using the uniform impor- 
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tance wavefunction resulted in larger statistical errors; 
unfortunately there is no easy way to construct a phase- 
separated trial function in this case. This is another ex- 
ample of the great care that must be taken in choosing 
an appropriate trial function in projector methodsEi In 
Table I we compare energies calculated using the SSE 
and CPMC methods for 32 site systems with U — 1 and 
V = 8. For the rest the results in this paper, we have 
used the SSE and IR methods exclusively, using SSE for 
energies when possible and IR for structure factors and 
susceptibilities. 



n 


IR 


SSE 


CPMC 


0.3125 


-0.461(2) 


-0.4644(2) 


-0.4656(2) 


0.3750 


-0.451(2) 


-0.4484(2) 


-0.4506(2) 


0.4375 


-0.382(3) 


-0.3825(1) 


-0.3829(2) 


0.5000 


-0.293(3) 


-0.3000(2) 


-0.3002(3) 


0.5625 


-0.211(3) 


-0.2105(2) 


-0.2108(4) 


0.6250 


-0.113(3) 


-0.1179(2) 


-0.1168(6) 


0.6875 


-0.025(4) 


-0.0261(2) 


-0.0200(6) 


0.7500 


-0.070(3) 


0.0629(3) 


0.072(1) 


0.8125 


0.150(2) 


0.1510(2) 


0.1640(7) 


0.8750 


0.239(2) 


0.2406(3) 


0.245(1) 



Table I: Comparison of QMC energies per site for V = 8 
and U = 1. Statistical errors in last digit are indicated 
in parentheses. The CPMC results used a free-electron 
trial function and were scaled for At — > from At = 0.1 
and Ar = 0.05. 



III. PHASE SEPARATION 
A. Phase Separation and Superconductivity 

We have already noted that an obvious motivation for 
mapping carefully the regions phase separation in the 
extended Hubbard model is that in a variety of strongly- 
correlated models, superconductivity has been found, or 
argued to be present, in close proximity to phase separa- 
tion boundary. Such proximity is intuitively reasonable, 
since both superconductivity and phase separation re- 
sult from effectively attractive interactions. The case of 
the EHM for V < provides an illustration. In this re- 
gion the model phase separates into a low-density phase 
and a high-density phase; the high-density phase con- 
sists either of adjacent doubly occupied sites or of adja- 
cent single electrons, depending on the value of [/J§E2l 
Near these phase-separated regions at negative V are 
well-established regions of singlet and triplet supercon- 
ducting fluctuations. 

However, the possible proximity of superconductivity 
to phase separation is also problematic, for if one applies 
the Luttinger liquid formalism and numerical techniques 
to calculate K p for relatively small systems in a region 
that is in fact phase separated in the thermodynamic 
limit, one will obtain incorrect results. Hence, before 
a mapping to Luttinger liquid parameters can be used, 



the phase separation boundaries of a model should be 
accurately determined. To quantify this point, we here 
note that in one recent study of the EHM at quarter 
filling)] the phase separation boundary for V 3> \U\ was 
determined using a cluster Gutzwiller approximation to 
find the vanishing of the inverse compressibility. The 
smallest V (for any of the values of U studied) for which 
phase separation was found was V ~ 14 (for U ~ — 2.5L 
An exact diagonalization study gave comparable results. □ 
Our QMC studies of larger systems reveal that phase 
separation occurs already at V ~ 8 at quarter filling, in 
close proximity to where these previous studies indicated 
that dominant superconducting fluctuations first appear. 

B. Strong— Coupling calculations in the V — > oo limit 

An .effective model, first derived by Penc and Mila 
(PM)J3 conveniently captures the exact behavior of the 
EHM for V — ► oo. For infinite V, any existing pairs can- 
not be broken up (or moved), while single electrons can- 
not occupy neighboring sites and hence behave as spinless 
fermions. The effective model thus consists of immobile 
pairs and single spinless electrons. For parameter val- 
ues where pairs and unpaired electrons coexist, the min- 
imum energy corresponds to having the pairs and the 
spinless fermions separated into two distinct regions, so 
as to minimize the kinetic energy of the spinless fermions 
by allowing them to move in the largest possible region. 
The pairs are separated by one lattice spacing. The en- 
ergy for a given number of pairs and single electrons can 
then readily be calculated analytically in the thermody- 
namic limit as the sum of U times the number of pairs 
plus the energy of a spinless fermion chain: 

mnU 2 , , . n 

E N = — £-(l-n)sm 7r (1 - m) . 2 

2 7r 1 — n 

Here m is the fraction of fermions forming pairs in the 
ground state, E/N is the energy per lattice site. The 
energy can be minimized with respect to m to determine 
the ground state. For U < —4, one can show that all elec- 
trons are paired for all fillings, giving a spin gap and the 
boundary of the Luther-Emery region for V — » oo. For 
U > —4, the ground state contains only unpaired elec- 
trons for fillings less than a critical filling, above which 
the pair phase starts to grow. At half filling, the sys- 
tem contains only pairs. Equation (^J) thus provides an 
exactly solvable model exhibiting phase separation. 

While PM focused exclusively on the quarter-filled 
case, their energy expression (||) can be used to obtain 
the phase-separated region for all fillings n. If the system 
is phase separated with the high- and low-density phases 
at densities 71hd and tild, respectively, then phase sep- 
aration occurs for the total (average) densities n in the 
range uld < n < njjD- In the thermodynamic limit the 
ground state energy must be linear as a function of n 
in this regime, reflecting the fact that adding particles 
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to the system only changes the relative sizes of the two 
phases in a system with fixed particle number (canonical 
ensemble) . 
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FIG. 1. The exact energy per site for for V — + oo and (7 = 1 
(solid circles, left axis) calculated using Eq. (|^). The depen- 
dence on n is completely linear for n > 0.37, as indicated by 
the dashed line, reflecting phase separation. The open squares 
(right axis) denote the fraction m of fermions forming pairs 
inEq. §). 

Figure [l] shows the energy per site versus filling in the 
V — > oo limit. At low fillings the system contains only un- 
paired electrons. Above a critical filling uld, the energy 
becomes linear, reflecting phase separation as pairs are 
formed. The linear behavior persists all the way up to 
half filling (n = 1 = tt-hd), because the high density phase 
consists of pairs occupying every other site. The phase 
separation density can be easily calculated as a function 
of U, resulting in the phase diagram shown in Figure |2|. 



1.0 



0.8 



0.6 



0.4 



0.2 



0.0 







B 






2? 


LU 


PHASE SEPARATION 


Degene 


PHAS 

red) 








MIFORM 

(All Pai 






UNIFORM PHASE 


=> 


/ (Spinless Fermion) 



-6.0 



-4.0 



-2.0 



0.0 

u 



2.0 



4.0 



6.0 



FIG. 2. The region of phase separation (PS) (extending 
from U = —4 to U = +4) for infinite V. The other phases 
are discussed in the text. 

For U > 4, the phase-separated and uniform states be- 



come degenerate at V — > oo, which we label as non-phase 
separated because we expect any finite V will break the 
degeneracy and the resulting state will not be phase sep- 
arated. Mila and Zotos havepshown that this state in 
fact is uniform and insulating.Q For finite V, we expect 
a smaller phase-separated region, as pairs will be able 
to break up with a finite amount of energy. Thus we 
are able to confine our numerical investigation of phase 
separation to the range —4 < U < 4. 



C. Numerical Calculation of the PS Boundary 

Numerically, for finite V , the phase separation bound- 
ary can be determined using the energy as a function 
of the number of particles in the system, as discussed 
above. However, a completely linear behavior will not 
be observed in a finite system, because the boundaries 
between the high- and low-density phases then occupy a 
significant fraction of the lattice and raise the energy per 
site by an amount that depends on the size of the bound- 
ary regions (we use periodic systems and hence have two 
boundaries). This will cause the energy as a function of 
n to be concave, which is not possible in the thermody- 
namic limit. A line can be drawn which is tangent to the 
E(n) curve at two points, which then constitute estimates 
of the fillings tt-ld and rtHD of the high- and low-density 
phases. In the absence of phase boundaries (which be- 
come irrelevant to the energy per site in the thermody- 
namic limit), the line corresponds to E/N of a system 
separated into regions of densities n L rj and rhd- For 
small system sizes, this "Maxwell construction" can be 
expected to be more accurate than signals of phase sepa- 
ration based on, e.g., probability distributions of particle 
densities, since the ground state energy typically exhibits 
only small finite-size effects in cases where the thermody- 
namic limit state is uniform. It_was previously assumed 
that rijjD = 1 also for finite Vu As we will see, in fact 
"hd < 1 in some parameter regions. In any case, we will 
refer to tild as the phase separation density. 

As a complementary method of determining the phase 
separation boundary, as well as to help in characterizing 
the phases, we have also used a criterion based on the 
static charge structure factor S p (q): 



S p (q) 



(3) 



Continuity at q — requires S p (q — ► 0) = in a uniform 
system. In a phase-separated periodic system of size N 
in the canonical ensemble, S p (q) will have a maximum 
at the shortest non-zero wave- number q = qi = 2ir/N, 
with a divergence as N — > oo. In small systems close 
to the boundary of phase separation, this signal is, how- 
ever, ambiguous, since there will be a range of parameters 
for which it is not possible to determine accurately the 
S(q —> 0) behavior based on the behavior for q > q%. We 
shall see examples of this in later sections. 



5 



Returning to the Maxwell construction, we note that 
for it to be reliable one must be certain that the energy 
can be computed accurately. Although there are no a 
priori approximations in two of the QMC methods we 
use, near half filling at large V, we have noticed that the 
QMC simulations may converge poorly, apparently be- 
coming stuck in meta-stable states. This occurs because 
the ground state near half filling contains a significant 
fraction of on-site pairs that require considerable energy 
to break up or move, making it difficult for the simulation 
to explore the full phase space of the model. The sim- 
ulations are particularly hard when the ground state is 
phase separated. The system can then separate into sev- 
eral alternating domains of high-and low-density phases, 
instead of just one of each. The resulting large statistical 
errors in the energy close to half filling cause problems in 
the determination of the tangent line and thus the point 
of phase separation. Fortunately, exactly at n = 1, the 
system is always uniform and the ground state energy 
can be calculated very accurately using exact diagonal- 
ization, since the finite size effects are very small at this 
filling (much smaller than for n < 1). This can be easily 
understood from perturbation theory around the static 
(t = 0) ground state, which is non-degenerate at n = 1 
(and only there) and consists of alternating doubly occu- 
pied and empty sites (i.e., a period- two CDW). A simple 
second-order perturbation calculation gives an energy per 
site of 
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E/N=-- — 

1 2 W-U 



(4) 



This thermodynamic expression compares remarkably 
well with exact diagonalizations on small systems; for ex- 
ample with U = 1 and V = 10, it gives E/N = 0.431034, 
while the exact results for N = 6 and N — 12 are 
E/N = 0.431092 and E/N — 0.431096. Unfortunately, 
the smallness of the finite-size effects holds only exactly 
at half filling (since the particles are no longer localized 
once the system is doped away from half filling) . 

For n < 0.7, the QMC simulations converge well even 
for phase-separated ground states (for V/t < 10), and 
we have been able to calculate E/N to within absolute 
statistical errors of ~ 10~ 4 — 10~ 3 using the SSE method 
(energies for small systems agree with ED results). Fig- 
ure H shows our Maxwell construction for N — 64, U = 1 
and V = 8. Also shown is a plot of AE, the difference 
between the E/N values and the best fit tangent line. 
Using this method, we find that the system in this case 
phase separates at n slightly above quarter filling. Data 
for N — 32 give the same result, indicating that this is 
indeed close to the phase separation point in the ther- 
modynamic limit. Thus, for a quarter filled system, we 
can conclude that Vps ~ 8.0 for U = 1 (since for V = 8, 
phase separation occurs at n only slightly above 1/2, and 
the critical filling decreases with increasing V). Previous 
studies at this fillingp'Q found a Vps ~ 14 - 18 at U = 1. 
Our phase separation boundary is hence at significantly 
lower V than reported and is very close to the K p = 1 



curve obtained using exact diagonalizationBEl This is a 
clear indication of possible problems with the prediction 
of a novel superconducting region for V 3> \U\. 
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FIG. 3. Maxwell construction for V = 8 and U = 1. The 
solid circles are energies calculated for N = 64. The dashed 
line goes through the n = 1 energy, and is tangent to the 
energy curve at the phase separation filling (density of the 
low-density phase) tild ~ 0.53. The inset shows the difference 
between the data points and the tangent line. 

A definite disadvantage of using the Maxwell construc- 
tion method for determining phase separation bound- 
aries is that the need to calculate very precise energies 
causes this approach to be extremely time consuming. 
Accordingly, as already noted, we employed also a sec- 
ond method of detecting phase separation, which relies 
on the behavior of the static charge structure factor de- 
fined in Eqn. (||) for q — » 0. Figure ^ shows S p (q) for 
several densities at U = —1, 0, 1, and V = 8. For U = 1 
we obtain phase separation at a density n « 0.59 from 
S p (q), compared to n w 0.53 from the Maxwell construc- 
tion. Note, however, that at U = 1 and n = 9/16, S p (q) 
shows no phase separation, whereas the Maxwell con- 
struction implies the system is already phase separated. 
S p (q) gives a slightly smaller phase-separated region than 
the Maxwell construction because a clear signal of an in- 
crease as q — > qi for a relatively small lattice requires that 
the system already be well inside the region of phase sep- 
aration. 

Examining the region of phase separation for U < 
reveals several interesting properties: whereas at U = 1 
the system stays phase separated for all densities n < 1 
above some critical density tild, for U = we see clear 
signs that phase separation disappears at high densities. 
This can be seen in the structure factors in Figure ^, 
which go smoothly to zero for U = —1,0 and n > 0.9. 
The Maxwell construction is also consistent with phase 
separation to a state with a high-density phase at filling 
less than n = 1, although the relatively large statistical 
errors in the energies close to half filling make it difficult 
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to obtain an independent estimate of bhd ■ However, in- 
vestigating the charge structure factor for several system 
sizes, it is clear that the uniform phase for n < 1 is sta- 
ble. The peak position of the charge susceptibility is also 
consistent with a high-density phase different from the 
[7 = 1 case (see next section). We conclude that there is 
a "re-entrant" uniform state for U ~ 0. In addition, from 
simulations carried out at larger V (V = 10), we see the 
phase-separated region moving further into the U < 
side of the phase diagram and the re-entrant uniform 
state remains. As we will discuss in the next Section, the 
re-entrant state has strong 2kp CDW fluctuations. 

We believe that the the re-entrance of the uniform state 
for U < 0, but not for U > 0, can be qualitatively under- 
stood in the following way: When V is large and n — ► 1, 
pairs are formed to avoid the nearest-neighbor repulsions. 
For U > 0, the pairs have positive energy. 
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FIG. 4. The static charge structure factor S p (q) for small 
momenta q/n at various fillings for V = 8 and (7 = —1,0, 1. 
Filled circles are for N = 32, open squares for TV = 64 and 
triangles for N = 128. As explained in Section IV, the straight 
lines all have slope = ir, corresponding to the q — + behavior 
if the LL exponent K p — 1. 



phase to be minimized to one lattice spacing (the kinetic 
energy of a pair is very small due to the V barrier to hop- 
ping, and the motion of the pairs can be be neglected). 
The combined effects of avoiding the repulsive on-site U 
and derealization of mobile unpaired electrons always 
favor phase separation as n — > 1 if V is large. 

On the other hand, for U < 0, the pairs have nega- 
tive energy and maximizing their number is favorable. 
For n close to 1 a uniform state with mobile pairs can 
then be stabilized if V is not too large. However, as the 
density is lowered the number of negative-energy pairs 
decreases, and if \U\ is small a phase separated state can 
then again be lower in energy due to the lower kinetic en- 
ergy of dominantly unpaired electrons in the low-density 
phase. Increasing \U\ (U < 0), maximizing the number 
of mobile pairs becomes more favorable, and the size of 
the region of phase separation decreases, in agreement 
with our results. 

For U > 1 at V = 8, we find that the phase separa- 
tion boundary moves towards higher fillings, becoming 
hard to discern for U ~ 4, as expected from the infinite 
V solution. For C7 < — 1, based on 32-, 64- and limited 
128-site data, phase separation does occur for V = 8, 
but is much harder to detect than at U = 1 and [7 = 0. 
Increasing V to 10, we cannot obtain as accurate results 
as for V = 8, but we do find clear signs that the phase 
separation region moves to lower values of n, as expected 
if there is a continuous evolution of the phase separation 
boundary to the infinite V solution discussed in the pre- 
vious section. We also find that phase separation does 
not occur at all for V = 4, and for V = 6 appears to be 
present only in a small region for U > 0. 

Finally, we stress the need to use more than one sys- 
n=1 1/16 tern size in applying the approach based on S p (q), as 
otherwise the analysis may give ambiguous or erroneous 
results for the phase separation boundary. For instance, 
based on the 32- and 64-site results shown in Fig. |[ there 
would appear to be no phase separation at U = — 1 and 
the slope of S p (q) versus q as q — > exceeds tt for some 
fillings. This would imply (see next section) that K p > 1 
and therefore dominant superconducting correlations in 
this region. However, based on data for larger (N = 128) 
lattices, it is likely that the system is phase separated at 
these fillings (see n — 9/16 for both U = — 1 and U = 
in Figure 0) . 



n=9/16 



n=15/16 



IV. CORRELATIONS AND FLUCTUATIONS 



A. Methods to calculate LL parameters 



Via phase separation, the total energy is minimized by 
balancing the potential energy of the pairs (which in- 
creases with the number of pairs) and the kinetic en- 
ergy of a low-density phase with dominantly unpaired 
electrons. The pressure of the unpaired phase forces the 
distance between the pairs in the immobile high-density 



In this section we describe results of calculations for 
a variety of correlation functions of the ground state of 
the EHM in non-phase-separated regions using various 
different estimators for the Luttinger liquid correlation 
exponents K p and K a . For models (like the EHM) with 



spin-rotation symmetry, the exponent K a = 



except 
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in the Luther-Emery phase, in which case only K p has 
meaning, since the system is spin-gapped and not a LL 
in the spin sector. For K a — 1, the asymptotic equal- 
time correlation functions in a LL have simple power law 
forms governed by K p : 

AWcdwM ~ r^ 1+K ^ cos(2fc F r), (5) 
^SS,TsM~r-( 1+1 /^>. (6) 

Hence, when K p < 1, SDW/CDW correlations domi- 
nate at large r, while for K p > 1, superconducting corre- 
lations dominate. At weak coupling, the renormalization 
group "g-ology" procecLurp_can be used to determine the 
exponents K p and K a P 4 tl 35 i Away from weak coupling, 
one must use the Luttinger liquid equations that relate 
K p to other physical observables of the .model. In our 
present notation, the relations we use areuJ 



K p = Dp/2v p , 
irv p K/2, 



K n 



K p = (itKDp/4) 1 / 2 , 



(7a) 
(7b) 
(7c) 



where D p is the Drude weight of the optical conductiv- 
ity (the charge stiffness), v p the velocity of the charge 
excitations, and k the compressibility. The use of three 
different relations for K p is important to verify the valid- 
ity the calculation. In particular, the three relations will 
give inconsistent results if the system is not a LL, or if 
finite-size or other systematic errors are present. 

The Drude weight can be G&Lc»lated from QMC simu- 
lations as the p — > oo limit oEJ'EZl 



D p = 7t[(-K/N) -A c (i27r/P)}, 



(8) 



where K is the kinetic energy, and K c (ioj m ) is the Mat- 
subara frequency-dependent charge current-current cor- 
relation function: 



1 f p 

A p (ioj m ) = — / dre % 
^ Jo 



'<?W(o)>. 



(9) 



The standard method used in exact diagonalization cal- 
culations of the compressibility k employs a finite differ- 
ence approximation involvings-ground state energies for 
different numbers of particles:ll3 explicitly, 



■ E(N, N e 



2E(N,N e )} 



(10) 

where E(N,N e ) is the ground state energy of N e elec- 
trons on an JV-site lattice. The compressibility may be 
also computed from the q — ► limit of the static charge 
susceptibility; 



where 



(«) = ij2 ei9U ' l) f" drfaWmiO)), (12) 
i,l 



where fij(r) is the charge at site j and imaginary time r. 
This definition avoids the errors due to discretization in 
the finite-difference definition (|l0|) , and we will therefore 
primarily use Eq. ( |ll|) here. Finally, the velocity v a as- 
sociated with the gapless charge (a = pi or spin (a = a) 
mode may be extracted from the raticO 



W a {q) = 2S a (q)/ X a(q), 



(13) 



where S a (q) is the static structure factor given by Eq. 
(3). W a (q) gives an upper bound to the energy of the 
lowest excitation of momentum q and becomes the ex- 
act excitation energy as q — > in a LL, so that v a can 
be directly extracted from the g-dependence for small q. 
Hence, the quantities needed for all three estimates of K p 
defined in Eqs. (0) can be computed directly from QMC 
data. Examining Eqs. (|7b|) , pi]), and (|l3|), one can see 
that K a is also given directly by the slope of the static 
structure factor as q — > 0: 



Kp. a = —Sp. a (q -> 0). 
irq 



(14) 



This relation may also be obtained directly from the 
calculation of the charsepcharge or spin-spin correlation 
functions in LL theory :ll3 



(UaOnar) 



K a cos(2kpr) 



(7rr) s 



,1+K a 



(15) 



« = xdq -> 0), 



(11) 



The Fourier transform of the non-oscillating term of 
equation ( |l5| ) leads to the expression (|lj) for K p and K a . 
As the structure factor is usually much easier to calcu- 
late numeric ally than the compressibility or the Drude 
weight, Eq. ([14|) is quite useful for calculating Luttinger 
liquid exponents. In particular, because S a (q) only de- 
pends on equal-time correlations in imaginary time, for 
most QMC methods the structure factors converge much 
faster than D p or the susceptibilities and hence may pro- 
vide a better estimate for K p . But the important caveat 
applies that the consistency among all three relations 
must be checked, which requires more detailed calcu- 
lations. Moreover, one additional consistency check is 
available in regions where one expects no spin gap. There 
it is required that K a — 1, which can be verified (or dis- 
proved!) from the q — ► limit of S a (q). Of course, the 
limit q — ► is impossible to attain strictly in numeri- 
cal simulations of finite systems. Since Eq. ( |l3| ) is an 
upper bound of the lowest excitation energy at momen- 
tum q, we expect that values of K p and K a calculated 
from equation ( p^ ) to be in general slightly larger than 
their true values. Effects of non-linearity in the true low- 
est excitation energy should be smaller for typical values 
of the smallest q accessible (i.e., we expect effect of the 
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broadening of the mode to be larger than effects of non- 
linearity) . 

One additional complication to using the slope of the 
structure factor to determine K p is that for very small 
fillings, 2kF may be close in momentum to the smallest 
q = q±. A strong broad CDW peak at this momentum 
may then affect the behavior all the way down to qi . The 
value obtained for K p will then still be an upper bound, 
but may be much larger than the true value. 

Since from Eq. (fry) one sees that K p is proportional 
to the compressibility, k, which (naively) should diverge 
at the phase separation boundary, calculations of K p us- 



ing Eq. (7b) or equivalently Eq. ( J14| ) near a region of 
phase separation typically show a strongly increasing K p . 
This increase of K p has often been interpreted (in a vari- 
ety of one-dimensional models) as evidence-fpr . ^.region 
of dominant superconducting correlationslj&Ej'OO. A 
more careful analysis begins by noting that in an infinite 
system, k jumps dis continuously from a finite value to in- 
finity at the phase separation boundary if the transition 
is first order, which is normally the case. However, as 
observed above ia-jSec. II, and also recently by Hellberg 
and Manousakis,E3 n in a finite system can diverge only 
inside the region where the infinite system is phase sepa- 
rated, due to the energy cost of the interface between the 
two phases. Hence, a result K p > 1 based on a diverging 
k obtained from small systems may be misleading, since 
it is possible that the system is already phase separated 
in the thermodynamic limit. 

As already noted, all these considerations raise legiti- 
mate concerns about the reliability of the previous cal- 
culations of the phase separation boundary in the EHM 
and hence also lead to concerns about whether K p ex- 
ceeds, unit y in the uniform phase. Indeed, a recent QMC 
studyliao of a two-band ID Hubbard model related to 
the EHM (a two-band "Cu-O" chain) showed that a first- 
order phase separation transition preempts superconduc- 
tivity in all but a small part oLthe phase diagram, in con- 
trast to a previous ED studyJlJ which indicated that su- 
perconducting fluctuations always dominate in the neigh- 
borhood of phase separation. 



B. QMC results for K p in the EHM 

Using QMC techniques to evaluate K p accurately via 
all three Eqs. (Q) is extremely time consuming, as D p in 
particular is difficult to calculate in the V U param- 
eter region. Accordingly, in this section we first present 
results using all three relations (Q) for only one value of U, 
and then use these detailed results to "benchmark" our 
alternate method of getting K p directly from Eq. (14). 
We choose U — 1 for these calculations, and find ex- 



128-site chains, we have performed Lanczos ED calcula- 
tions for 16 sites, in order to investigate systematically 
the finite-size effects. We have also checked the QMC 
simulations against the ED data for this system size. In 
our ED calculations, we use Eq. ( |To| ) to define the com- 
pressibility k and to extract the charge velocity from the 
lowest charge excitation energy. 

We begin with a comparison (for U — 1 and varying V) 
of the various low-energy parameters used in calculating 
K p . Figure || compares the Drude weight, the charge 
velocity, and the compressibility as computed for n = 1/2 
(quarter filling). Up to V w 2 the ED and QMC results 
agree almost perfectly. For larger V the deviations are 
due to the larger finite-size effects in the ED data. For 
V < 6, the finite-size errors in the ED results seem to be 
greatest for v p ; hence we expect that Eq. ([ft]), which is 
the only Luttinger liquid relation not involving v p , will 
give the best estimate for K p in small systems. For large 
V, the slower increase of the compressibility as computed 
by ED is likely primarily due to discretization errors in 
Eq. [l^, which become large in regions where the energy 
curvature changes rapidly as a function of N e (as is the 
case close to phase separation). 




cellent agreement between the two methods, which then 
justifies our use of the static structure factor method for 
other values of U. 

In addition to the QMC simulations for 32-, 64-, and 



FIG. 5. Comparison of the Drude weight (D p ), the charge 
velocity (v p ), and the compressibility (ft) for U — 1 as com- 
puted by ED with N = 16 (solid curves) and QMC with 
N = 64 (triangles). Where not shown, the QMC error bars 
are approximately the size of the symbols (somewhat larger 
for V = 7, 8). The dashed lines provided guides to the eye for 
the QMC data. 

Figure [?] shows results for K p as a function of V cal- 
culated for U = 1 by Lanczos diagonalization and QMC 
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simulations. The data points based on Eq. (7b) are equiv- 
alent to using only the structure factor in determining K p 
according to Eq. (14). The error bars of the other meth- 
ods are not shown for clarity; in general they are at least 
twice as large as the errors using the structure factor 
slope method. Further, in Figure 0, the K p values for 
V = 8 calculated using Eqs. ( |7a| ) and ( |7c| ) have been left 
off, as their statistical errors are too large for this large 
V. Given the equivalence of using Eq. ([7l]) to the calcu- 
lation of the slope of the structure factor at q = 0, it is 
clear that the slope method provides a good estimate for 
K p , as it agrees with the other methods within error bars 
but shows soi&ller statistical fluctuations. Previous K p 
calculationsQtO using ED methods had seen large differ- 
ences among the different relations (0). As seen in Fig. |], 
these differences are greatly reduced using larger lattices. 
Indeed, our QMC results on 64-site systems show that 
the three relations give equivalent results, to within er- 
ror bars. Most importantly, all the K p values obtained 
for V < 8 and U = 1 are less than one. The Maxwell 
construction, shown in Fig. |^, indicates that phase sepa- 
ration occurs at quarter filling for V only slightly larger 
than V = 8. It is therefore clear that there is no extended 
region where K p > 1 before phase separation, although 
there is a definite increase as the phase separation bound- 
ary is approached. We can of course not strictly rule out 
the existence of an extremely narrow region where K p 
exceeds one. However, this seems unlikely in view of our 
results for the spin susceptibility at V = 8 and 10 shown 
in Figure ^. In both cases there is a clear peak at q — 2kp, 
which would not be expected if K p exceeds unity. Since 
our Maxwell construction indicates that the V = 10 
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FIG. 6. x<y for U = 1 and n = 0.5 (quarter-filling) at V = 8 
(solid circles) and V — 10 (open circles). 

system should be phase separated in the thermodynamic 



limit, we can conclude that the dominant SDW fluctua- 
tions persist at the phase separation boundary and there 
is no region with K p > 1. 
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FIG. 7. K p vs. V for U = 1 calculated using all Eqs. (Q). 
ED results (curves) for 16 sites are compared with QMC re- 
sults for 64 sites (symbols). Error bars are shown only for 
Eq. (7b) for clarity; error bars for Eqs. (jT^) and ([7c]) are larger. 
Phase separation for this filling and U occurs at V ~ 8 

For other fillings and U we have concentrated on 
V = 8. For 1 < U < 4 we find similar behavior in K p 
(based on the slope of S p (q) for 32 and 64 site systems) 
to the U = 1 case: K p increases as the density increases, 
but does not exceed one before phase separation. For 
—3 < U < 0, our determination of the phase separation 
boundary is not as accurate, and in some cases calcu- 
lations for for 32- and 64-site systems give K p > 1 in 
regions where there are no clear signals of phase separa- 
tion. However, in cases where we have also carried out 
calculations for N = 128, the general trend seems to favor 
phase separation over a stable uniform state with K p > 1 
(recall the data for U = — 1 and U — at n — 9/16 in 
Fig. |). 

The persistence of finite-size effects in calculations of 
Vps and K p close to phase separation for systems as large 
as N = 64 sites emphasizes the importance of studying 
and understanding finite-size effects in calculations of K p 
and related quantities. In Sec. V we will discuss the com- 
plete phase diagram of the model and comment further 
on the behavior of K p as V — > Vps- 



C. Unusual charge correlations in the EHM 

Luttinger liquid theory predicts structure in the charge 
or spin response only for multiples of 2kp, which is a 
consequence of the low-energy effective model being lin- 
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earized about ±kp. Peaks are typically seen only at 
q = 2/cp and/or q = 4fcp (only for the charge if there 
is a spin gap). The 2kp peaks should diverge as T -^Jl 
if K p < 1, and 4fcp peaks are divergent for small if p .t£l 
Divergences in finite systems can of course be seen only 
as N — > oo. 

In studying the V ^> U region of the EHM, we have 
found strong charge response peaks also at wavevectors 
that are not multiples of 2kp\ we shall use the term 
"anomalous peaks" to refer to these unusual charge cor- 
relations. Importantly, however, these anomalous peaks 
appear to be associated with gapped modes. As the tem- 
perature decreases and N increases they do not diverge. 
Hence we believe the system is still a LL in this region 
for sufficiently low energies. We shall comment further 
below on the interpretation of these anomalous peaks. 

Figure || shows the charge susceptibility for several fill- 
ings for V = 8 and U = —1, 0, 1. Starting with the U = 1 
data, we see that at small fillings the susceptibility is 
dominated by a large peak at q — Akp. This corresponds 
to a system with almost no pairs, where the particles be- 
have essentially as spinless fermions. For higher fillings 
the large peak shifts from the value of 4fcp to a slightly 
lower momentum — which for later purposes we call 4fc F ff 
— whose value depends on the filling as well as the other 
model parameters. As one moves further into the phase- 
separated 
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FIG. 8. Charge susceptibility for V = 8, U = -1,0, 1. All 
data are for 64 site systems. The vertical solid and dashed 
lines indicate q = 2fcp and q = 4fep, respectively. Note that 
for n > 1/2, 4&f = 7r(2 — 2n) in the reduced zone scheme. 
Lines between data points are guides to the eye. 

region (above n ks 5/8) a strong peak develops at q = w, 
independent of n (and hence kp). This corresponds to 
the wave-vector of the high-density phase in the phase- 



separated region, which is a CDW state in which double 
occupancies alternate with empty sites. 

Turning to the cases for U — — 1 and U — 0, the be- 
havior is quite similar for low and intermediate densities. 
As in the U = 1 case, at low fillings one sees a large 
4fcp peak and a much weaker (but still observable in the 
U = data) 2k p peak. For intermediate fillings near the 
phase separation boundary, the main peak is no longer 
at 4&f but rather at 4fc F as for U = 1. This is shown 
clearly in the U = data for n = 7/16, which is outside 
of the phase-separated region (n < tild)- At large fill- 
ings (n > tihd), instead of the peak at q = tt seen in the 
U = 1 data, the U = — 1 data shows a strong peak at 
2kp , reflecting the re-entrance (as a function of n) of the 
homogeneous phase, which is characterized by 2k-p CDW 
fluctuations. A weaker 4/cp peak can also be seen at this 
filling. For U = 0, n = 7/8 is very close to our estimated 
re-entrance point (tihd ~ 0.88). 

Figure ^ shows the location of the dominant peak in 
the charge susceptibility as a function of filling for V = 8 
and U — — 1,0, 1. For small fillings, all the systems are 
dominated by the Akp charge response. As the filling is 
increased, the anomalous 4A: F ff peak dominates the re- 
sponse for these system sizes. For still higher fillings the 
normal 2k? peak dominates. As already noted, we see a 
peak also at the densities where the anomalous peak 
is present (see Figs. || and 0), and we expect it to 
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FIG. 9. Momentum value of the largest peak in the charge 
susceptibility Xp(q) f° r a 64-site system as a function of band 
filling n for V = 8, U = 1,0,-1. The peak position was 
determined by fitting a 2 nd degree polynomial to points near 
the peak. The solid and dotted lines indicate q — 2kp and 
q = 4&f respectively. 

diverge as N ~ + oo and T ~ * 0, reflecting asymptotic LL 
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behavior. Note the constant behavior of the peak posi- 
tion for U — at 0.65 <■ n < 0.88; this indicates that this 
peak originates from the high-density phase of the phase- 
separated system, which apart from growing in size as n 
increases remains unchanged in the phase separation re- 
gion. This behavior is seen less clearly for U = — 1, where 
the region of phase separation is apparently very narrow. 
Close to half filling a uniform state again stabilizes for 
U = —1,0, whereas the U = 1 system remains phase 
separated all the way up to half filling. 

Since the anomalous charge response is not at a 
wavevector possible within the Luttinger liquid formal- 
ism, it is important to determine whether the correspond- 
ing 4fcp ff peak diverges with decreasing temperature and 
thus has thermodynamic relevance in the low-energy sec- 
tor. Figure ^ shows the momentum dependence of the 
charge susceptibility at three different temperatures for a 
parameter set where the anomalous charge fluctuations 
are strong. The inset shows the peak value versus in- 
verse temperature for two different system sizes. We see 
that as the temperature is lowered, the anomalous peak 
does not diverge (it in fact has a maximum at a finite 
temperature), and there is almost no dependence on N. 
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FIG. 10. Behavior of non-LL "anomalous" 4fcp charge 
peak in Xp(q) for V = 8, (7 = —1, n = 3/8 vs. temper- 
ature. The symbols connected by dashed lines in the main 
figure are for 64-site systems. The inset shows the peak value 
vs. (3 for 32-site (circles) and 64-site (squares) systems. Note 
also the presence of a normal 2&f peaks which grows with 
increasing /3 

Nonetheless, this non-LL charge peak is still relatively 
low in energy. Using equation ( |T3| ) we are able to estimate 
an upper bound of E ~ 0.75i above the ground state 
at V = 8, U = -1 and n = 3/8. In Figure [[(], one 
clearly sees that the normal 2kp peak indeed grows with 



decreasing temperature, as expected in a Luttinger liquid 
with K p < 1, although the amplitude is relatively weak 
even at /3 = 256. 

We believe that the anomalous charge response is due 
to on-site pairs that are sufficiently long-lived (due to the 
difficulty of breaking them up via high-energy intermedi- 
ate states) to form a CDW together with the single parti- 
cles of the system. Specifically, in the infinite V effective 
model of Penc and Mila,Q there are two kinds of parti- 
cles; spinless fermions and bosons with charge two. The 
number of doubly occupied sites in the simulation of the 
EHM corresponds to the number of bosons, and hence 
we can calculate the average total number of particles 
(fermions plus bosons) within the effective model. We 
find that the momentum of the anomalous charge peak 
is consistent with spacing uniformly a number of particles 
equal to this total the number of particles. Hence, there 
is a tendency to form a CDW consisting of a mixture 
of fermions and bosons which repel each other. Never- 
theless, the non-divergence of the anomalous peak, the 
presence of a 2k f peak, and the consistency among the 
Luttinger liquid relations indicate that the asymptotic 
low-energy properties are still those of spin- 1/2 fermions 
forming a LL. 



D. Spin gap in the EHM 

Our results for the spin susceptibility of the extended 
Hubbard model in the V \U\ region show either a 
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FIG. 11. X"(q) an d Sail) comparing spin-gapped to 
non-gapped cases. Left panels are for V = 8, U = — 1, and 
n = 0.375, where no gap is present. Right panels are for same 
U and V but n — 0.75, where the system is gapped. Solid 
curves connecting data points are guides to the eye. The dot- 
ted lines indicate slope n for 5 CT (g) vs. q. 

normal Luttinger liquid 2k p peak or the presence of a 
spin gap. In general, we find that the spin response is 
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much noisier and harder to converge numerically than 
the charge response, as expected in a parameter region 
dominated by charge correlations. To determine whether 
the system is spin gapped, we examine the spin suscep- 
tibility Xa-io) m the limit q — > 0. If Xu(q — ► 0) = 0, 
the system is spin gapped. The presence of a spin gap 
can also be inferred from the structure factor S a (q): If 
there is no spin gap, then K a = 1, which translated into 
a slope 7r for the structure factor versus q (as discussed 
in Sec. IVA). In a spin-gapped system S a {q) should ap- 
proach zero faster than linearly and in particular should 
fall below the line with slope tt (instead of approaching 
this line from above, as expected if there is no spin gap). 
Again, one must use care in phase-separated regions, as 
the x<?(<l) wm show a mixture of responses from both 
phases. Figure |Tl| shows examples of X<t{q) anc l S a (q) for 
parameters with and without a spin gap. The two differ- 
ent behaviors discussed above can be clearly discerned. 
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FIG. 12. Summary of spin gap boundary. The solid curve 
is the g-ology result for V = 1. Hollow symbols connected 
by dotted lines are results of ED calculations and finite-size 
scaling for the values of V indicated. Solid symbols connected 
by lines are QMC results for 32 or 64 sites. 

When interpreting the complete ground state phase 
diagram, it is useful to consider the development of the 
spin gap boundary as V increases from V = 0. For the 
case V = 0, there is a spin gap (and dominate singlet su- 
perconducting fluctuations) for all negative U. For weak 
coupling, g-ology results may be used, and they predict 
a spin gap for g\ = U + 2Vcos(7m) < 0. For larger U 
and V we have used both QMC and Lanczos ED. In the 
ED calculations, the gap was computed from the ground 
state energies of the system, and the system with one spin 
flipped: A s = E{N,n^,n{) - E(N,n 1 + l,n; - 1). For 
each filling, two or three system sizes and finite-size scal- 
ing of the gap values versus 1/N were used. The QMC 
estimate for the spin gap boundary was extracted from 
structure factor results such as those shown in Fig. [ll| . 



We consider a system to be gapped in cases where the 
S a (q) curve falls below the line with slope ir as q — > 0. 
The QMC and ED results are combined in Figure |lj. 
One can see that (i) for weak interactions the ED results 
are quite close to the g-ology predictions; and (ii) impor- 
tantly, as V increases, the spin-gap boundary line moves 
towards the infinite V result of U = —4. 
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FIG. 13. The spin gap region as a function of U and V for 
7i = 1/4, n = 1/3, n = 1/2, and n — 2/3. The dotted lines in- 
dicate the phase boundaries from g-ology, and the solid circles 
are the results of our exact diagonalizations. The spin gap is 
present to the left of the boundaries shown. The approxima- 
tive rfi>fease separation boundaries for negative U calculated in 
Refsat£l are also shown. 

For purposes of comparison to previouSrjesults on the 
spin gap in the parameter region V > Ojj'u we present 
ED results graphed in the (U, V) plane for various fill- 
ings in Figure [l^. In this figure we also show the g- 
ology predictions for the dominant fluctuations, as well 
as the phase separation boundary for the V < regime, 
which was calculated previously^ Our spin gap boundary 
for n = 1/2 agrees, closely with previous numerical work 
presented by PM.Q However, our boundary for n = 2/1 
is significantly different from that obtained previously,!!! 
with our results placing the spin gap boundary further 
towards the negative-U side of the phase diagram. 



V. DISCUSSION AND CONCLUSION 

In the preceding sections, we have explored numeri- 
cally the V ^> \U\ ~ 1 parameter region of the ID ex- 
tended Hubbard model for a wide range of band fillings. 
Our numerical approaches included three different forms 
of QMC simulations on systems of up to 128 sites and 
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Lanczos exact diagonalizations and finite-size scaling for 
systems of up to 16 sites. In addition, we relied on some 
analytic results for V — ► oo. 
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FIG. 14. Phase diagram of EHM in the region -4 < U < 4. 
Circles denote exact diagonalization data. The curves for 
V = 0, V = 1, and V = 2 are from g-ology. At V = 8, 
filled boxes are phase separation points we have confirmed 
[via a Maxwell construction or the q — > behavior of S p (q) 
using QMC]. Open boxes are also most likely correspond to 
phase separation, but our results here are less certain. The 
full phase-separated region is inside a curve enclosing all these 
points. The curves shown in panel 5 are only schematic. The 
curves in panel 6 are the exact solution for V — > oo 



A complete quantitative mapping of the phase diagram 
even for the restricted region V \U\ ~ 1 would require 
enormous numerical work to determine precisely all phase 
boundaries for many values of the parameters. Although 
we have not carried out such a program in the present 
study, we believe that we have developed an accurate 
qualitative (with quantitative results for specific param- 
eters) picture of the phase diagram in a parameter region 
of the EHM about which previously little has been un- 
derstood. The considerable complexity we have found in 
the model serves to illustrate the great care that must be 
taken in interpreting numerical data, in particular those 
used as signals of superconductivity and phase separation 
in exact diagonalization results for small lattices. 

Our most important conclusion is that phase separa- 
tion extends to much lower values of V than previously 
reported. As a result, the Luttinger liquid exponent K p 
does not exceed one before phase separation, and hence 
the ground state is not dominated by superconducting 
fluctuations. This resolves the difficulties in interpre- 
tations of previous exact diagonalization results J3 which 
indicated an extended region of gapless spin excitations 
and K p > 1, which taken together would correspond to 
dominant triplet pairing correlations. Since the naive pic- 
ture of superconductivity in the V ^> \U\ region involves 
singlet on-site pairs, one would rather have expected a 
spin gap to accompany K p > 1. Our results explain this 
puzzling result as simply due to difficulties in detecting 
phase separation in small systems. It should be noted 
that our results for K p in the uniform phase agree quite 
well with the previous estimatesErQ up to our calculated 
phase separation boundary. 

To summarize our results for the EHM phase diagram 
for the V ^> \U\ region, it is most useful to examine how 
the features (phase separation, spin gap, possible super- 
conductivity, etc., evolve together from the more familiar 
V ~ region toward V = oo. Further, since exact results 
exist in both the V = and V — oo limits, following the 
evolution of the various boundaries among the phases 
from the two known solutions is very helpful in under- 
standing the global evolution of the phase diagram. Fig- 
ure |l4| combines our numerical data, the exact limiting 
cases, and qualitative considerations to summarize our 
results. The numbered comments below correspond to 
the panels of figure 14, labeled from top (1) to bottom 
(6): 

1. (V = 0) The EHM here reduces to the standard 
Hubbard model and is spin gapped with dominant 
superconducting fluctuations for all U < 0. 

2. (V = 1) As V increases from zero, the single 
(vertical) boundary that for V=0 separates the 
superconducting/spin-gapped phase from the Lut- 
tinger liquid/non-spin-gapped splits into two sepa- 
rate boundaries, producing a central spin-gapped 
but non-superconducting region, which extends 
into both positive and negative U. A region of 
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singlet superconductivity is also present. Results 
of g-ology agree well with exact diagonalization at 
this weak coupling. 

3. (V = 2) For this intermediate value of V, 
the "standard" superconducting region has disap- 
peared from our parameter region, moving to larger 
negative U values. The spin gap region also recedes 
further towards negative U . 

4. (V = 8) At large V, a region of phase separation 
has entered. Based on the movement of the phase- 
separated region towards negative U as V increases, 
phase separation probably first enters near the end 
of the spin gap line, i.e. near n = 1 and between 
U = and U = 4. For U < 0, the spin gap bound- 
ary remains largely unchanged from weak coupling. 
At positive U, the high-density phase consists of 
pairs on every other site. For negative U, the high- 
density phase expands slightly, as the dominant 
wavevector is no longer at q — n. In addition, 
the region of phase separation exists over a lim- 
ited region of n and vanishes near half filling, with 
the homogeneous phase re-entering. For small fill- 
ings, 4fcp charge correlations dominate, while near 
half filling 2kp charge correlations dominate. At 
intermediate fillings we find an anomalous, non- 
Luttinger liquid peak that is non-divergent with 
increasing system size or decreasing temperature. 
This peak appears to reflect CDW fluctuations in- 
volving long-lived pairs that repel each other as well 
as the unpaired particles. Our data, which clearly 
show phase separation at this V, show no convinc- 
ing signs of dominant superconducting correlations 
in the uniform phase, leading us to conclude that 
any region of superconductivity for V 3> \U\ is very 
small indeed, or, more likely, does not exist. 

5. (V 3> U) At still larger V, we expect the region 
of phase separation to grow towards the point U — 
—4, n = 0. We show a schematic phase diagram in 
this case. 

6. (V — oo) The point connecting the spin gap and 
phase separation boundaries has moved to (U, n) = 
(—4, 0). The spin gap boundary is now the vertical 
line from n = 0ton = lat?7 = —4. The phase 
separation region consists of either effective spinless 
fermions or pairs, and the high density phase has 
pairs on every other site. 

One may argue that our numerical data cannot exclude 
the existence of an extremely narrow strip of supercon- 
ductivity preceding phase separation. However, in the 
non-spin-gapped region this seems unlikely in view of the 
obvious difficulty in accounting for dominant triplet su- 
perconducting fluctuations in the V 3> \U\ region. The 
exact V — ► oo solution also provides a counter- argument: 
below the phase separation boundary (n < jt-ld), the 



system is mapped onto spinless fermions, and the dom- 
inant correlations are then clearly not superconducting 
in the neighborhood of the phase separation boundary. 
For large but finite V, we find numerically that K p does 
not appear to exceed unity close to the phase separation 
boundary. We instead find strong 2&f spin fluctuations 
that do not vanish as the phase separation boundary is 
crossed. Hence we believe that the system is dominated 
by SDW fluctuations adjacent to the phase separation 
boundary in cases where there is no spin gap. On the 
other hand, for U < we find K p « 1 in a narrow re- 
gion for which we cannot definitely conclude that the 
system is phase separated. This ambiguous region coin- 
cides with our calculated spin gap boundary, and is essen- 
tially the region of the open boxes in panel 4 of Fig. [T^ . 
Hence, if superconducting fluctuations indeed dominate 
in this region, they would be of the singlet type, which 
is what one would expect. However, also in this case 
it appears more likely that K p does not exceed unity 
at the phase separation boundary. The V — > oo case 
again provides support for this result: On the left-hand 
side of the phase-separated region shown in panel 6 of 
Fig. |lj (i.e., U < —4), the system contains only pairs, 
and for large but finite V exhibits strong 2fcp CDW fluc- 
tuations. Dominant superconducting fluctuations do ap- 
pear in this region for small values of V but are replaced 
by CDW fluctuations for moderate V (see panels 2 and 3 
of Fig. |lj). It then appears unlikely that a small region 
of dominant superconducting fluctuations in the region 
U ~ (—3, —1) and n ~ (0.2,0.5) would re-appear as V is 
increased and then again vanish as V — > oo, but we do 
not have data that can definitely exclude this scenario. 

Although we believe our results establish that phase 
separation, rather than superconductivity, dominates the 
V 3> \U\ region of the EHM, there remain a number of 
interesting open questions about some of the "exotic" 
non-Luttinger liquid effects we have observed. In partic- 
ular, two such questions involve (i) the exact nature of 
the state corresponding to the anomalous 4fc^ peak in 
the charge susceptibility and (ii) a more detailed under- 
standing of the "re-entrant" phase separation behavior 
compared to the normal U > phase separation, where 
the high-density phase is at (the naively expected) half- 
filled density. We have presented qualitative interpre- 
tations of these "exotic" effects, but they may warrant 
further study, focusing on when they can appear, what 
simple effective model (if any) can be used to describe 
them, the relation between phase separation and the non- 
Luttinger liquid peak, and how this peak evolves with 
system size and temperature. One possible method to ac- 
cess this parameter region would be to study numerically 
an effective model including paired and single electrons 
similar to PM's V — * oo effective model. Such an effec- 
tive model would eliminate electrons on nearest-neighbor 
sites, which could be additionally taken into account per- 
turbatively, if necessary. We are currently exploring this 
approach. 
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